;;----------------------------------------------------------------
;; calculate the percentile of vertical profiles   
;;
;; kai zhang, max planck institute for meteorology 
;; 
;;----------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "callvar.ncl"

begin 


  fna = "BC_SP2_TC4_2007_WB57"      ;; ng/m3 
  fna = "BC_SP2_CR-AVE_2006_WB57"   ;; ng/m3 
  fna = "BC_SP2_ARCTAS-2008_DC8"    ;; ng/m3  
  fna = "BC_SP2_ARCTAS-2008_P3B"    ;; ng/m3 STP 
  fna = "BC_SP2_ARCPAC-2008_WP-3D"  ;; ng/m3 
  fna = "BC_SP2_ICE-L_2007.nc"      ;; ng/m3 
  fno = "profile_"+fna+".nc"  
 
  fna = fna + ".nc" 
  fla = addfile(fna,"r") 

  system("rm "+fno) 
  flo = addfile(fno,"c") 
 
  aa = fla->BC 
  bb = fla->BCP  

  d1 = dimsizes(aa)
  d2 = dimsizes(bb(0,:)) 

  pp = new((/5,d2/),"float") 


  do ik = 0,d2-1 

     ;;-----------------------------------------------------------------
     ;; remove missing values 
     ;;-----------------------------------------------------------------

     print("level    :  "+ik) 

     xx = bb(:,ik) 

     if(.not.(all(ismissing(xx)))) then 

     cc = vnmiss_1D(xx) 
     qsort(cc) 

     ;;-----------------------------------------------------------------
     ;; get 10% 25% 50% 75% 90% percentile index 
     ;;-----------------------------------------------------------------

     dm = dimsizes(cc) 

     x10 = round(.10*dm,3)-1
     x25 = round(.25*dm,3)-1     ; -1 to account for NCL indexing starting
     x50 = round(.50*dm,3)-1     ; -1 to account for NCL indexing starting
     x75 = round(.75*dm,3)-1     ; at 0            
     x90 = round(.90*dm,3)-1

     pp(0,ik) = (/ cc(x10) /) 
     pp(1,ik) = (/ cc(x25) /) 
     pp(2,ik) = (/ cc(x50) /) 
     pp(3,ik) = (/ cc(x75) /) 
     pp(4,ik) = (/ cc(x90) /) 

     delete(cc) 

     end if 
     
  end do 


  pp!0 = "pt" 
  pp&pt= ispan(1,5,1) ;; (/"10%","25%","50%","75%","90%"/)
  pp&pt@_1 ="10%"
  pp&pt@_2 ="25%"
  pp&pt@_3 ="50%"
  pp&pt@_4 ="75%"
  pp&pt@_5 ="90%"
  pp!1 = "lev"
  pp&lev = ispan(1,d2,1)  
  pp@long_name = "percentile" 
 
  flo->profile=pp 


  do ip = 0,4

     print("=== ip "+ip+" ===") 
     print("") 
     print("") 
     print("") 
     print(""+pp(ip,:)) 

  end do 



end


